Fitting GP to lightcurve 428_…in Stan

Notebook outlining the fitting of GP to thunderKAT lightcurve ID$ 428_…

1 Light Curve

  • The light curve has \(N =\) 21 observations over a range of 147.57 days.
  • Observations are evenly spread over the time range.
  • The shortest gap between observations is 5 days.
  • The longest gap between observations is 14.68 days.
  • The mean flux density is \(\bar{y} =\) 0.19 Jy.
  • The mean standard error in the observations is 0.000491 Jy.
  • The observational noise is very small relative to the brightness of the observations.

2 SE Basic Model

  • Zero constant mean function.
  • Squared Exponential kernel function.
  • Homoskedastic noise.
  • Wide prior on observational noise, uninformed by observational noise estimates.

\[y \sim \mathcal{N}(f(x), \sigma_\textrm{noise}^2)\]

\[f \sim \mathcal{GP}(\boldsymbol{0}, k(x, x'))\]

\[k(x,x') = \eta \exp\left\{ -\frac{1}{2}\frac{(x - x')^2}{\ell^2}\right\}\]

\[\ell \sim \mathrm{InvGamma}(5,5)\]

\[\eta \sim \mathcal{N}^+(0,1)\]

\[\sigma_\textrm{noise} \sim \mathcal{N}^+(0,1)\]

2.1 MCMC Results

Warning: 2 of 4000 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
$num_divergent
[1] 2 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 1.043763 1.221051 1.102958 1.218542
 variable     mean   median      sd     mad      q5      q95   rhat ess_bulk
    eta     0.0960   0.0421  0.1648  0.0386  0.0092   0.3781 1.0007     2447
    ell   162.1374 139.9640 90.2990 67.1896 66.0926 330.1677 1.0013     2901
    sigma   0.0063   0.0061  0.0012  0.0011  0.0047   0.0085 1.0000     3299
 ess_tail
     2595
     2530
     2310

2.2 MCMC Plots

2.3 Posterior Predictive Samples

The fitted model has a very long lengthscale, comparable to the length of the observational window. The estimated observational noise has a standard deviation more than an order of magnitude of that in the original data. The combination of these parameters has lead to a very smooth fit that passes through the middle of the observed data points rather than through any datapoints themselves.

2.4 PSD

3 SE Observational Errors Model

  • Zero constant mean function.
  • Squared exponential kernel function.
  • Heteroskedastic (Gaussian) noise.
  • Incorporate data on error in observations of each \(y_i\).

\[y_i \sim \mathcal{N}(f(x_i), \sigma_i^2)\]

\[f \sim \mathcal{GP}(\boldsymbol{0}, k(x, x'))\]

\[k(x,x') = \eta^2 \exp\left\{ -\frac{1}{2}\frac{(x - x')^2}{\ell^2}\right\}\]

\[\ell \sim \mathrm{InvGamma}(5,5)\]

\[\eta \sim \mathcal{N}^+(0,1)\]

\[\sigma_i \sim \mathcal{N}^+(\textrm{stderr}(y_i), \mathrm{Var}(\textrm{stderr}(\boldsymbol{y})))\]

$num_divergent
[1] 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 0.9342605 1.0233746 0.9841971 1.0595353

3.1 MCMC Results

 variable     mean   median      sd     mad       q5      q95    rhat ess_bulk
 eta       0.02216  0.01994 0.01039 0.00778  0.01081  0.04093 1.00042     6217
 ell      11.22759 11.23455 0.51423 0.49400 10.35169 12.04060 1.00085     6382
 sigma[1]  0.00039  0.00039 0.00004 0.00004  0.00033  0.00046 1.00137     8109
 ess_tail
     2839
     2855
     2866

3.2 MCMC Plots

3.3 Posterior Predictive Samples

By including the observed observational errors for setting priors on the Gaussian noise of each observation, the fitted median passes through each of the observed points.

3.4 PSD

4 Matern 3/2 kernel

  • Matern 3/2 covariance kernel
  • Zero constant mean function

\[y_i \sim \mathcal{N}(f(x_i), \sigma_i^2)\]

\[f \sim \mathcal{GP}(\boldsymbol{0}, k(x, x'))\]

\[k(x,x') = \eta \left( 1 + \frac{\sqrt{3(x - x')^2}}{\ell}\right) \exp\left\{ -\frac{\sqrt{3(x - x')^2}}{\ell}\right\}\]

\[\ell \sim \mathrm{InvGamma}(5,5)\]

\[\eta \sim \mathcal{N}^+(0,1)\]

\[\sigma_i \sim \mathcal{N}^+(\textrm{stderr}(y_i), \mathrm{Var}(\textrm{stderr}(\boldsymbol{y})))\]

$num_divergent
[1] 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 1.023875 1.080972 1.061123 1.047926

4.1 MCMC Results

 variable     mean   median       sd      mad       q5       q95    rhat
 eta       0.04154  0.02482  0.06405  0.01692  0.00897   0.12613 1.00033
 ell      86.05902 79.21425 30.28509 22.88986 52.29066 141.91955 1.00006
 sigma[1]  0.00039  0.00039  0.00004  0.00004  0.00032   0.00046 1.00113
 ess_bulk ess_tail
     3175     2357
     3121     2219
     7029     3063

4.2 MCMC Plots

4.3 Posterior Predictive Samples

4.4 PSD

5 SE + Matern 3/2 additive kernel

  • Sum of squared exponential and Matern 3/2 kernels
  • single output scale (marginal variance) hyperparameter
  • zero constant mean function

\[y_i \sim \mathcal{N}(f(x_i), \sigma_i^2)\]

\[f \sim \mathcal{GP}(\boldsymbol{0}, k(x, x'))\]

\[k(x,x') = \eta \left[ \exp\left\{ -\frac{1}{2}\frac{(x - x')^2}{\ell_\mathrm{SE}^2}\right\} + \left( 1 + \frac{\sqrt{3(x - x')^2}}{\ell_\mathrm{M}}\right) \exp\left\{ -\frac{\sqrt{3(x - x')^2}}{\ell_\mathrm{M}}\right\} \right]\]

\[\ell_\mathrm{SE} \sim \mathrm{InvGamma}(5,5)\]

\[\ell_\mathrm{M} \sim \mathrm{InvGamma}(5,5)\]

\[\eta \sim \mathcal{N}^+(0,1)\]

\[\sigma_i \sim \mathcal{N}^+(\textrm{stderr}(y_i), \mathrm{Var}(\textrm{stderr}(\boldsymbol{y})))\]

5.1 MCMC Results

$num_divergent
[1] 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 1.0384751 1.0393340 0.9615068 0.9812617
 variable     mean   median       sd      mad       q5      q95    rhat
 eta       0.01194  0.00898  0.01151  0.00477  0.00394  0.02844 1.00093
 ell_SE   40.27136 35.43965 20.98407 13.58655 19.76977 76.49121 1.00026
 ell_M    55.55175 52.95730 15.53669 12.78883 35.85322 82.75450 0.99958
 sigma[1]  0.00039  0.00039  0.00004  0.00004  0.00032  0.00046 1.00001
 ess_bulk ess_tail
     3476     2008
     4438     2183
     3255     2159
     5398     2705

5.2 MCMC Plots

5.3 Posterior Predictive Samples

5.4 PSD

6 SE + Matern 3/2 (2 output scales) additive kernel

  • Sum of squared exponential and Matern 3/2 kernels
  • One output scale (marginal variance) hyperparameter for each kernel
  • zero constant mean function

\[y_i \sim \mathcal{N}(f(x_i), \sigma_i^2)\]

\[f \sim \mathcal{GP}(\boldsymbol{0}, k(x, x'))\]

\[k(x,x') = \eta_\textrm{SE}^2 \exp\left\{ -\frac{1}{2}\frac{(x - x')^2}{\ell_\mathrm{SE}^2}\right\} + \eta^2_\textrm{M}\left( 1 + \frac{\sqrt{3(x - x')^2}}{\ell_\mathrm{M}}\right) \exp\left\{ -\frac{\sqrt{3(x - x')^2}}{\ell_\mathrm{M}}\right\}\]

\[\ell_\mathrm{SE} \sim \mathrm{InvGamma}(5,5)\]

\[\ell_\mathrm{M} \sim \mathrm{InvGamma}(5,5)\]

\[\eta_\textrm{SE} \sim \mathcal{N}^+(0,1)\]

\[\eta_\textrm{M} \sim \mathcal{N}^+(0,1)\]

\[\sigma_i \sim \mathcal{N}^+(\textrm{stderr}(y_i), \mathrm{Var}(\textrm{stderr}(\boldsymbol{y})))\]

Warning: 5 of 8000 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
$num_divergent
[1] 3 0 1 1

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 1.0236362 0.9947824 1.0108117 0.9614624

6.1 MCMC Results

 variable     mean   median       sd      mad      q5      q95   rhat ess_bulk
 eta_SE     0.0000   0.0000   0.0000   0.0000  0.0000   0.0001 1.0006     6224
 eta_M      0.0357   0.0183   0.0692   0.0136  0.0057   0.1101 0.9998     8500
 ell_SE     1.2993   1.0819   0.8226   0.4905  0.5507   2.7425 1.0003     9555
 ell_M    317.2728 258.4510 228.5529 151.3097 97.7173 744.3742 1.0003     6149
 sigma[1]   0.0004   0.0004   0.0000   0.0000  0.0003   0.0005 1.0002    13074
 ess_tail
     3197
     5063
     4418
     5179
     5494

6.2 MCMC Plots

6.3 Posterior Predictive Samples

6.4 PSD

7 SE x Matern 3/2 multiplicative kernel

\[y_i \sim \mathcal{N}(f(x_i), \sigma_i^2)\]

\[f \sim \mathcal{GP}(\boldsymbol{0}, k(x, x'))\]

\[k(x,x') = \eta \exp\left\{ -\frac{1}{2}\frac{(x - x')^2}{\ell_\mathrm{SE}^2}\right\}\left( 1 + \frac{\sqrt{3(x - x')^2}}{\ell_\mathrm{M}}\right) \exp\left\{ -\frac{\sqrt{3(x - x')^2}}{\ell_\mathrm{M}}\right\}\]

\[\ell_\mathrm{SE} \sim \mathrm{InvGamma}(5,5)\]

\[\ell_\mathrm{M} \sim \mathrm{InvGamma}(5,5)\]

\[\eta \sim \mathcal{N}^+(0,1)\]

\[\sigma_i \sim \mathcal{N}^+(\textrm{stderr}(y_i), \mathrm{Var}(\textrm{stderr}(\boldsymbol{y})))\]

Warning: 72 of 8000 (1.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
$num_divergent
[1]  0  0 72  0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 0.9895645 0.9546933 0.8873692 1.0126744

7.1 MCMC Results

  variable     mean   median       sd      mad      q5       q95    rhat
 eta        0.03333  0.02631  0.04699  0.01695 0.00984   0.07239 1.12222
 ell_SE    41.47464 41.01155 33.60317 25.98761 0.73175  97.00669 1.52568
 ell_M     57.84199 65.94690 38.30645 26.86249 0.73907 110.92935 1.52388
 sigma[1]   0.00039  0.00039  0.00004  0.00004 0.00033   0.00046 1.00042
 f_star[1]  0.18392  0.18392  0.00040  0.00039 0.18326   0.18458 1.00020
 ess_bulk ess_tail
       22     2881
        7       27
        7       28
     7949     5107
     7611     7813

7.2 MCMC Plots

7.3 Posterior Predictive Samples

7.4 PSD

8 SE + Periodic

\[y_i \sim \mathcal{N}(f(x_i), \sigma_i^2)\]

\[f \sim \mathcal{GP}(\boldsymbol{0}, k(x, x'))\]

\[k(x,x') = \eta \left[ \exp\left\{ -\frac{1}{2}\frac{(x - x')^2}{\ell_\mathrm{SE}^2}\right\} + \exp\left\{ -\frac{2 \sin^2\left( \pi\frac{\sqrt{(x - x')^2}}{T}\right)}{\ell_\mathrm{P}^2}\right\} \right]\]

\[\ell_\mathrm{P} \sim \mathrm{InvGamma}(5,5)\]

\[\ell_\mathrm{SE} \sim \mathrm{InvGamma}(5,5)\]

\[\eta \sim \mathcal{N}^+(0,1)\]

\[T \sim \mathcal{U}[\textrm{minimum gap in x}, \textrm{range of x}]\]

\[\sigma_i \sim \mathcal{N}^+(\textrm{stderr}(y_i), \mathrm{Var}(\textrm{stderr}(\boldsymbol{y})))\]

8.1 MCMC Results

Warning: 26 of 8000 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
$num_divergent
[1]  1 12  3 10

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 0.9677225 1.0269724 0.9582524 1.0597445
 variable    mean  median      sd     mad      q5      q95   rhat ess_bulk
   eta     0.0036  0.0032  0.0018  0.0013  0.0017   0.0068 1.0026     2771
   ell_SE  9.9347  9.9336  0.7520  0.6500  8.7242  11.1683 1.0043      581
   ell_P   2.7280  2.4112  1.3876  0.9705  1.2502   5.2363 1.0046     1307
   T      92.5493 95.0402 34.0438 40.6806 36.5269 142.2803 1.0023     1305
 ess_tail
     2963
      174
      257
      513

8.2 MCMC Plots

8.3 Posterior Predictive Samples

8.4 PSD

9 SE + Matern 3/2 + Periodic

\[y_i \sim \mathcal{N}(f(x_i), \sigma_i^2)\]

\[f \sim \mathcal{GP}(\boldsymbol{0}, k(x, x'))\]

\[k(x,x') = \eta \left[ \exp\left\{ -\frac{1}{2}\frac{(x - x')^2}{\ell_\mathrm{SE}^2}\right\} + \left( 1 + \frac{\sqrt{3(x - x')^2}}{\ell_\mathrm{M}}\right) \exp\left\{ -\frac{\sqrt{3(x - x')^2}}{\ell_\mathrm{M}}\right\} + \exp\left\{ -\frac{2 \sin^2\left( \pi\frac{\sqrt{(x - x')^2}}{T}\right)}{\ell_\mathrm{P}^2}\right\} \right]\]

\[\ell_\mathrm{P} \sim \mathrm{InvGamma}(5,5)\]

\[\ell_\mathrm{SE} \sim \mathrm{InvGamma}(5,5)\]

\[\ell_\mathrm{M} \sim \mathrm{InvGamma}(5,5)\]

\[\eta \sim \mathcal{N}^+(0,1)\]

\[T \sim \mathcal{U}[\textrm{minimum gap in x}, \textrm{range of x}]\]

\[\sigma_i \sim \mathcal{N}^+(\textrm{stderr}(y_i), \mathrm{Var}(\textrm{stderr}(\boldsymbol{y})))\]

9.1 MCMC Results

Warning: 7 of 8000 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
$num_divergent
[1] 7 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 0.8965104 0.9769791 0.9785691 0.9511891
 variable     mean   median      sd     mad      q5      q95   rhat ess_bulk
   eta      0.0041   0.0033  0.0030  0.0016  0.0016   0.0089 1.0150      291
   ell_SE  25.8580  23.3266 15.3052  8.9439  1.7806  50.8321 1.0716       37
   ell_M   35.7989  36.0480 13.1031  9.5932  2.0152  55.2239 1.0716       38
   ell_P    2.5098   2.2009  1.2929  0.8982  1.1765   4.8095 1.0035     3120
   T      108.3781 114.5820 29.7569 30.1309 51.7676 144.4434 1.0137      300
 ess_tail
     2822
       11
       11
     3895
      104

9.2 MCMC Plots

9.3 Posterior Predictive Samples

9.4 PSD